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Abstract 

This paper presents rigorous forward error bounds for linear conic optimization 
problems. The error bounds are formulated in a quite general framework; the un- 
derlying vector spaces are not required to be finite-dimensional, and the convex cones 
defining the partial ordering are not required to be polyhedral. In the case of linear pro- 
gramming, second order cone programming, and semidefinite programming specialized 
formulas are deduced yielding guaranteed accuracy. All computed bounds are com- 
pletely rigorous because all rounding errors due to floating point arithmetic are taken 
into account. Numerical results, applications and software for linear and semidefinite 
programming problems are described. 



1 Introduction 



In this paper forward error bounds for the optimal value of linear conic optimization problems 
as well as certificates of feasibility and infeasibility are presented, including the discussion of 
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rounding-off errors and details of implementation. These rigorous bounds aim to prove how 
accurate the approximate results computed by any conic solver are. The underlying vector 
spaces are in general infinite-dimensional, that is the bounds are developed in the framework 
of functional analysis. 

Forward and backward error analysis together with a detailed discussion of rounding-off 
errors and condition numbers for matrix problems were first described in the outstanding 
papers published sixty years ago by von Neumann and Goldstine [32] and Turing [45]. Turing 
writes: 

Error estimates can be of two kinds. We may wish to know how accurate a 
certain result is, and be willing to do some additional computation to find out. 
A different kind of estimate is required if we are planning calculations and wish to 
know whether a given method will lead to accurate results. In the former case we 
do not care what quantities the error is expressed in terms of, provided they are 
reasonably easily computed. With these estimates we wish to be absolutely sure 
that the error is within the range stated, but at the same time not to state a range 
which is very much larger than necessary. With the second type of estimate, the 
error is preferably expressed in terms of quantities whose meaning is sufficiently 
familiar that the general run of values involved may at least be guessed at. 

Particularly, forward error bounds for the inverse of a matrix including a discussion of the 
effects of rounding-off errors are presented there. Today one would speak in this context of 
verified or rigorous error bounds, and thus these two papers can be viewed as the pioneering 
work in the field of verification methods, a part of numerical analysis. Forward error bounds 
are propagated in interval arithmetic; see the textbooks Alefeld and Herzberger pQ, Moore 
|22j . and Neumaier [27], [28]. But also in other areas the interest in rigorous forward error 
bounds is growing. Parlett [31], for example, remarks in relation to the numerical accuracy 
of eigenvalue problems: 

For some of us, however, it has taken nearly 40 years to realize that backward 
stability is not enough. 

Also Trefethen writes in [H] about the future of Numerical Analysis 

I expect that most of the numerical computer programs of 2050 will be 99% 
intelligent and just 1% actual "algorithm" if such a distinction makes sense. 
Hardly anyone will know how they work, but they will be extraordinarily powerful 
and reliable, and will often deliver results of guaranteed accuracy. 
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Linear conic optimization refers to problems with a linear objective function and linear 
constraints where the variables are restricted to a cone. In general these problems are 
non-smooth. Linear programming, quadratically convex programming, second order cone 
programming and semidefinite programming are special cases. Since each convex problem 
can be described equivalently as a linear conic problem, the latter provides a universal 
form of convex programming (see Nemirovski [21], [25] ). Thus not surprisingly, a large 
variety of applications of conic programming are known from areas like system and control 
theory, combinatorial optimization, signal processing and communications, machine learning, 
quantum chemistry, and many others. For an elaborate bibliography the reader is referred 
to Wolkowicz [Tf] . 

Nesterov and Nemirovski [25] have shown that self-concordant barriers apply to many conic 
problems yielding polynomial time interior point methods. Renegar j3U] has investigated 
the sensitivity of infinite-dimensional conic optimization problems, and in [37] he analyzed 
interior point methods. He introduced a condition number for conic optimization, which is 
a generalization of the condition numbers defined by von Neumann, Goldstine, and Turing. 
This condition number is the scale-invariant reciprocal of the smallest data perturbation that 
will render the perturbed data instance either primal or dual infeasible. It is used in sensi- 
tivity analysis and moreover can be viewed as a problem instance size of conic optimization 
problems yielding important results in complexity theory. A problem is called ill-posed if 
this condition number is infinite, that is the distance to primal or dual infeasibility is zero. 

One of Renegar's main results is that the sensitivity of the optimal solutions and the optimal 
value can be bounded by the condition, and especially he proved that the bounds for the 
optimal value depend cubically on the inverses of the relative distances to primal and dual 
infeasibility. Renegar shows that this bound cannot be improved in general. For an ill-posed 
problem this result means that there exist arbitrarily small perturbed data instances such 
that the difference between the optimal value of the original problem and the perturbed 
problem is arbitrarily large, but the optimality conditions for the perturbed problem almost 
coincide with the optimality conditions for the original problem. Since conic solvers are 
terminated when the optimality conditions are satisfied approximately it cannot be distin- 
guished between the optimal values of the original and the perturbed problem in the case of 
ill-conditioned or ill-posed problems. A consequence is that the noise introduced by floating 
point arithmetic may occasionally yield to wrong termination and nonsensical computational 
results. 

In this paper we show how certain weak boundedness qualifications on e-optimal solutions 
can be used to compute rigorous forward error bounds for the exact optimal value, also for 
ill-conditioned or even for ill-posed problems. Such qualifications and even more restric- 
tive assumptions, like certain smoothness properties, are customary when solving ill-posed 
problems with regularization methods. It need not to be assumed that Slater's constraint 
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qualifications are fulfilled. The rigorous error bounds provide more safety for conic optimiza- 
tion problems, and they provide rigorous results in branch-and-bound algorithms for global 
and combinatorial optimization problems. Another application are computer-assisted proofs 
where it is mandatory to control all rounding-off errors (see for example Neumaier [30J and 
Rump [H]). It should be made clear that we do not investigate regularization methods. In 
this paper we assume that approximations computed by some conic solver (with or without 
regularization) are given, and these approximations are then used for computing the error 
bounds. It is of particular importance that the computation of the error bounds can be done 
outside the code of any imaginable solver as a reliable postprocessing routine, providing a 
correct output for the given input. Especially, we show for some combinatorial problems 
how branch-and-bound algorithms can be made safe, even if ill-posed relaxations are used. 
Numerical results for some ill-posed and ill-conditioned problems are included. 

Ill-conditioned and ill-posed problems are not rare in practice, they occur even in linear 
programming. Ordonez and Freund [33] stated that 71% of the lp-instances in the NETLIB 
Linear Programming Library [26] are ill-posed. This library contains many industrial prob- 
lems. Recently Freund, Ordonez and Toh 2006 [8] have shown that 32 out of 85 problems of 
the SDPLIB are ill-posed. 

The presented results in this paper formalize a viewpoint which apparently has not been 
made in conic programming. They can be viewed as an extension of methods for linear 
programming ( [T3] and Neumaier and Shcherbina [31]), and for smooth convex program- 
ming (see [T3]) to ill-conditioned and ill-posed non-smooth problems using the framework of 
functional analysis. 

The paper is organized as follows. After introducing some notation and basic definitions 
in Section 2, we consider in the next section conic optimization problems. Then in Section 
4 verified lower and upper bounds of the optimal value in the infinite-dimensional case are 
presented, and applied to finite-dimensional linear programming problems. Sections 5 and 
6 are devoted to error bounds for second order cone and semidefinite programming, respec- 
tively. Then in Section 7 we investigate conic optimization problems with block structured 
variables. In Section 8 verified certificates of infeasibility are presented, and in the following 
section we will focus on some applications in combinatorial optimization. Section 10 con- 
tains numerical results for the NETLIB Linear Programming Library (obtained by the C++ 
software package Lurupa p2]) and for the SDPLIB benchmark problems (obtained by the 
MATLAB software package VSDP [IS]). Finally, some conclusions are given. 
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2 Notation and Preliminaries 



Let X be a real vector space equipped with a norm ||.||, and let DC C X be a convex cone, i.e. 

DC + DC C DC, «3( C X for a G R+, (1) 

where R + denotes the set of nonnegative real numbers. A convex cone DC induces a partial 
ordering 

x<y<=>y-x<E%, (2) 

which is a transitive and reflexive binary relation on X compatible with addition and scalar 
multiplication: 

x < y, u < v, a G R+ =>- x + u < y + v and ax < ay. (3) 

Conversely, each partial ordering determines a convex cone, namely the positive cone 

DC := {x G X : x > 0}. (4) 

A vector space X equipped with a partial ordering is called a partially ordered vector space. 
A partial ordering is called antisymmetric, if 

x < y, y < x =>- x = y. 

It can be proved that antisymmetric partial orderings correspond to pointed cones, i.e. 

dc n (-DC) = {o}. 

If not explicitely mentioned we do not assume that the partial ordering is antisymmetric. 
Given a partial ordering the set 

[x, x] := {x G X : x < x < x} = (x + DC) n (x - X) (5) 

is called an interval. For a subset M of a partially ordered vector space X a vector x is called a 
lower bound of M, if x < m for all m G M, and in this case we write x < M. The lower bound 
x is called infimum of M if every other lower bound y of M satisfies y < x. Analogously, upper 
bounds and supremum are defined. X is said to be a vector lattice for the partial ordering < 
if for all x, y G X the supremum sup{x, y} and the infimum infja;, y} exists and is contained 
in X, respectively. In a vector lattice the operations x + := sup{x,0}, x~ := inf{x, 0} and 
|x| := sup{x, — x} are defined, and the properties |x| = x + — x~ , x = x + + x~ , \x\ = iff 
x = 0, | Ax | = |A| |x| for real A, and \x + y\ < \x\ + \y\ are satisfied. 

Let X* denote the dual space of X, that is the space of continuous linear functionals endowed 
with the operator norm. The set DC* of all positive linear functionals, i.e. 

DC* = {y G X* : (y, x) := y(x) > for all x G DC}, (6) 
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is a convex cone in X* denning a partial ordering in the dual space. 

The basic properties and relations for vector lattices as well as examples can be found in 
Birkhoff [3] and Bourbaki [5]; see also Peressini [35| and Schaefer [33]. We use the same 
notation |.| and < for all norms and partial orderings. It will always be clear from the 
context which norm and which cone is referred to. Hence, if x G X then x > means i6 3C, 
and if y G X* then y > denotes y G X*. Observe that we do not write y* for a continuous 
linear functional in X* because from the position in (y, x) the meaning is clear, and we can 
omit the star. This notion is closely related to Hilbert spaces and the Theorem of Riesz 
which states that the continuous linear functions can be represented by the inner product 
(y, x) where y is a vector in the Hilbert space. 

In the following some illustrative and well-known examples of normed vector lattices are 
shown. The real finite dimensional space X = R n equipped with the Euclidean inner product 
and the Euclidean norm |.| can be ordered by the convex cone 

X := R™ = {x G R n : Xi> for i = 1, . . . , n}. (7) 

This cone is self-dual (i.e. X = X*) and implies the lattice operations 

xf = max{0,Xj}, x~ = min{0,Xj}, \xi\ = xf — x^ (8) 

for i — 1, . . . , n. This vector lattice is used in linear programming (LP). 

In second order cone programming (SOCP) the same normed space X = R n is equipped with 
the partial ordering defined by the convex ice-cream or Lorenz cone 

% ■= j x = ^ : ^J e R n : x n > ||x : ||J , (9) 

where x. := (xi, . . . ,x n ^i) T . This cone is also self-dual and further properties are described 
in Section 5. 

In semidefinite programming (SDP) the real linear space X is R^+i)/ 2 , which is identified 
with the set of real symmetric n x n matrices X. The inner product of X, Y is defined 
by (X,Y) := traceX T y = E^-X^-Y^-, and the induced norm ||X|| := (traceX T X)s is the 
Frobenius norm. 

The space X = R, n (™+ 1 )/ 2 i s a Hilbert space, thus self-dual, and it is equipped with the 
self-dual cone of positive semidefinite matrices 

% := S 1 ™ = {X G X : X is positive semidefinite}. (10) 

Using the eigenvalue decomposition X = Q T AQ of a real symmetric matrix it follows that 

X- = Q T A-Q, X + = Q T A + Q, \X\ = Q T \A\Q, (11) 
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where A , A + , and |A| denote the diagonal matrices with nonpositive, nonnegative, and 
modulus of the eigenvalues of X on the diagonal, respectively. 

Occasionally, it is useful to represent symmetric matrices X as column vectors x by using 
the svec operator: 

x := svec(X) := (X 11 ; V2X 21 ; . . . ; s/2X nl ; X 22 ; V2X 32 ; . . . ; X nn ). (12) 

Here we follow the convention of MATLAB and use ";" for adjoining vectors in a column. 
Then the inner product between symmetric matrices X and Y is the usual inner product, 
that is 

(X,Y)=x T y. (13) 

We also use the notation x G % and x < y if the corresponding symmetric matrices X and 
Y such that x = svec (AT) and y = svec(F) have these properties. 

For any compact Hausdorff space Q, the vector space X := C(Q) of real- valued functions is 
a normed vector lattice with norm 

||/|| c( o) :=sup{|/(:r)|} 

and ordering cone 

% ■= {f G C(Q) : f(x) > for all x G Q}. 

Finally, we mention L P (Q), the vector space of Lebesgue-integrable functions / : Q — > R, 
where Q C R n , 1 < p < oo. This space is equipped with the norm 

\\f\\ P := (/ \f(x)\ P dx)K 
h 

and can be partially ordered by the cone 

% := {/ G L P (Q) : f(x) > almost everywhere on Q}. 

This yields a normed vector lattice, which is of interest in the case of Volterra and Fredholm 
type equations. 

3 Conic Optimization Problems 

We study rigorous error bounds for the conic optimization problem in standard form 

minimize (c, x) s.t. Ax = b, x G X, (14) 

where X is a real normed vector space, OC C X is a convex cone, c G X*, ^ is a real normed 
vector space, A denotes a continuous linear operator from X to y, and 6 G ^. With f p 



7 



we denote the primal optimal value, where f p := +00 if the problem is infeasible. Many 
interesting examples of optimization problems can be formulated in this framework. In the 
following some familiar facts are described. 

The Lagrangian function of problem (|14|) has the form 

L(x,y) := (c,x) + (y,b- Ax), (15) 

where y £ The optimization problem 

inf sup L(x, y) (16) 

xGOC y& y* 

is equivalent to ( |T4l) . Indeed, if b — Ax = then (y, b — Ax) = for each y G V*, and the 
supremum of L(x,y) is equal to (c, x). In the case where b — Ax 7^ there is some y with 
(y, b — Ax) > 0, and hence the supremum is infinite. 

Obviously, the Lagrangian satisfies L(x,y) = (y,b) + (— A*y + c,x) where A* is the adjoint 
operator. By exchanging in ( fT6l) infimum and supremum we obtain the dual problem 

sup inf L(x, y) (17) 

with optimal value fd- Since exchanging inf and sup always produces a lower bound, weak 
duality holds, that is fd < f p . Because inf L(x, y) = —00 whenever —A*y + c ^ X* the dual 
problem can be written equivalently in the form 

maximize (y, b) s.t. - A*y + c G %*, y G T- (18) 

We set fd '■= —00, if the dual problem is infeasible. 

Let x be primal feasible, and let y be dual feasible, then 

(c, x) = (c, x) + (y,b- Ax) = {-A*y + c, x) + (j/, b) > (y, 6), (19) 

and hence equality holds iff the complementarity condition 

(-A*y + c,x) = (20) 

are fulfilled. This condition means that the feasible pair x, y is a saddle point of the La- 
grangian. Moreover, it follows that there is no duality gap between the primal and the dual 
problem, and both problems have optimal solutions if and only if there exists a primal and 
a dual feasible solution fulfilling the complementarity conditions. In other cases where such 
primal dual feasible pairs do not exist strong duality may be not fulfilled. 

Duality theory is central to the study of optimization. First, algorithms are frequently based 
on duality (like primal-dual interior point methods), secondly they enable to check whether 
a given feasible point is optimal, and thirdly it allows to perform a sensitivity analysis. For 
more results on duality theory in the infinite-dimensional case see for example Renegar 
Rockafcllar 1 38 1 , and the literature cited there. 
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4 Lower and Upper Bounds for the Optimal Value 



This section is elementary but important for understanding both, the basic ideas behind 
rigorous forward error bounds and implementations. It turns out that for computing these 
error bounds only approximate primal and dual solutions x and y are required. Further 
assumptions about the accuracy of the approximations are not necessary; they need to be 
neither primal nor dual feasible. If the accuracy is poor, however, then the error bounds 
cause overestimation. 

The cones % and %* create partial orderings for the vector spaces X and X*, respectively. 
We assume in this paper that for subsets of these partially ordered vector spaces there exist 
lower and upper bounds. If the existence of the infimum or the supremum is necessary we 
mention this explicitely. Note that all vector lattices satisfy this assumption. 

We start with a simple result concerning bounds for linear functionals. 

Lemma 4.1 Assume that x, x G %, x < x, and let d,d~ G X* with dT < {d, 0}. Then 

(d,x) > (d~,x) and (d~,x) < 0. (21) 

Proof. Since < d - d~ and x > it is 

< (d — d~ , x) = (d, x) — (dT, x). 

Hence (d, x) > (d~ ,x). Since — dT > and x — x > the linearity of the functional — d~ 
yields 

< (— d~ , x — x) = {d~ , x) — {d~ , x) , 
which immediately implies (J2TT) . □ 
The dual version of this lemma is: 

Lemma 4.2 Assume that x,x~ G X with xT < {x, 0} , and let d,d G %* with d < d. Then 

(d,x)>(d,xT) and (d,x~) < 0. (22) 



Proof. Since < x — x~ and d > it is 

< (d, x — x~) = (d, x) — (d,x~). 
Hence (d,x) > (d,x~). Since — x~ > and d — d > it follows that 

< (d — d, —x~) = {d,xT) — (d,x~) 
which implies ( l2"2"j) . □ 
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Lemma 4.3 Assume that X and X* are normed vector lattices. Assume that for x,x G X 
and d,d G X* £/ie inequalities 

\x\ < x and \d\ < d (23) 

are satisfied. Then 

\(d,x)\ < (d,x) < \\d\\ \\x\\. 
Proof. Since \x\ = sup{x, — x} < x it follows that 

—X < X < X. 

Analogously we obtain 

-d < d < d. 
The inequalities x — x > and d — d > imply 

< (d — d, x — x) = (d,x) — (d, x) — (d, x) + (d, x) , 
and the inequalities x + x > and d + d > imply 

< (d + d, x + x) — (d, x) + (d, x) + (d, x) + (d, x) . 
Adding both inequalities yields 

< 2((d,x) + (d,x)), 

and therefore 

— (d, x) < (d, x). 

Because of the symmetry of x and — x in the definition \x\— sup{x, — x} it follows that 

(d,x) = —(d,—x) < (d,x). 

Hence 

\(d,x)\ < (d,x), 

and the last inequality follows from the definition of the norm of an operator. □ 

For bounding rigorously the optimal value, we claim boundedness qualifications, which are 
more suitable for our purpose than Slater's constraint qualifications. We assume that the 
conic optimization problem satisfies the following condition which we call primal boundedness 
qualification (PBQ): 

(i) Either the primal problem is infeasible, 

(ii) or f p is finite, and there is a simple bound x G % such that for every e > there exists 
a primal feasible solution x(e) satisfying x(e) < x and (c, x(e)) — f p < e 
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Observe that PBQ implies that the primal problem is bounded from below, but the existence 
of an optimal solution is not demanded, only simple bounds x for ^-optimal solutions are 
required. The following theorem provides a finite lower bound / of the primal optimal 
value. 



Theorem 4.1 Assume that PBQ holds. Let y £ ^* and let d := —A*y + c. Suppose further 
that dT < {d, 0}, then: 

(a) The primal optimal value is bounded from below by 

l>(y,b) + (dr,x)=:l p (24) 

(b) If dT = 0, then y is dual feasible and fa > / = (y,b), and if moreover y is optimal 
then Jd = f_ v - 

Proof, (a) If the primal problem is infeasible, then f p = +oo, and each finite value is a 
lower bound. Hence, assume that PBQ (ii) is satisfied with x := x(e) and e > 0. Then 

(c,x) = (d,x) + (A*y,x) 

= (y,b) + (A*y,x)-(y,b) + (d,x) 
= (y,b) + (y,Ax-b) + (d,x). 

Since x is primal feasible Ax — b = 0, and 

(c,x) = (y,b) + (d,x). 

Lemma 14.11 implies the inequality 

(c,x) > (y,b) + (gT,x). 

Because of PBQ (ii) 

f P > (c,x) - e> (y,b) + (d~ ,x) -e. 
For e — ► the assertion (a) follows. 

(b) If d~ — then d £ %*, implying that y is dual feasible, and the assertion follows. □ 

In particular, an approximate solution y which is close to optimality implies that d is close 
to %*. Hence, each lower bound d~ sufficiently close to d~ is almost zero, and it follows that 
f Pd (y, b) is reasonable; that is the overestimation is not very much larger than necessary. 
The lower bound uses the approximate optimal value (y, b) , and a correction is added which 
takes into account the violation of dual feasibility d~ evaluated at the upper bound x. 
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We illustrate the bound for linear programming problems in standard form 



minimize c T x s.t. Ax — b, x > 0. 



(25) 



This is the special case of the conic optimization problem where X = X* = R n , % = %* = R™ 
and y = = R m . It is well-known that for linear programming strong duality f p = fd ='■ f 
holds without any constraint qualifications. Hence, Theorem 14.11 yields immediately the 
lower bound 



of rounding errors for computing / by using directed rounding or interval arithmetic. The 
MATLAB toolbox INTLAB [12] provides the directed rounding modes, and the following 
short INTLAB program produces a rigorous lower bound: 

setround(-l) ; 

dlminus = min(0,A , *(-yt)+c) ; 
flow = b'*yt + dlminus ' *xup ; 
setround(O) ; 

If interval arithmetic is used, then the input data A, b, c may be intervals, and we obtain 
a lower bound for each instance within the interval data. Verified error bounds for general 
linear programming problems also with free variables can be found in [T3], and for formula 
(j2Bj) see Corollary 6.1 in [H]. 

To compute a rigorous upper bound of the optimal value we assume that the conic optimiza- 
tion problem satisfies the following condition, which we call the dual boundedness qualification 



(i) Either the dual problem is infeasible, 

(ii) or fa is finite, and there is a simple bound y such that for every e > there exists a 
dual feasible solution y(e) satisfying \y(e)\ < y and fa — (y(e), b) < e 

Theorem 4.2 Assume that DBQ holds. Let iel, and suppose further that 




n. It is straightforward to control all effects 



(26) 



(DBQ): 



\Ax — b\ < r, 

xT < {x, 0}, and 

d > —A*y + c for all dual feasible y with \y\ <y. 



(27) 
(28) 
(29) 



Then: 
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(a) The dual optimal value is bounded from above by 

f d <(c,x)-(d,x-) + (y,r)=:J d . (30) 

(b) If x~ = and r = 0, then x is primal feasible and f p < f d = (c, x), and if moreover x 
is optimal, then f p = f d . 

Proof, (a) If the dual problem is infeasible then = — oo, and each finite value is an 
upper bound. Hence, assume that DBQ (ii) is satisfied with y := y(e) and e > 0. Then 

(y,b) = -(y,Ax-b) + (y,Ax) 

= (c, x) + (y, Ax) - (c, x) - (y, Ax - b) 
= (c, x) — (c — A*y, x) — (y, Ax — b) . 

Since y is dual feasible and d > d := —A*y + c > 0, we can apply Lemmata 14.21 and 14.31 
which yield 

(y,b) < (c,x) - (d,x~) + (y, f). 

Because of DBQ (ii) 

fd < (y,b) + s < (c, x) - (d, x~) + (y, r) + e. 
For e^Owe obtain the upper bound (j3"0j) . 

The assertion (b) follows immediately, since xT = and r — 0. □ 

Observe that for finite dimensional y or in the case where ^ is a vector lattice the absolute 
value |.| is defined. If the absolute value is not available, then we replace the inequalities 
\y{z)\ < V, \Ax — b\ < f by ||y(e)|| < y and \\Ax — 6|| < r, respectively, and we obtain the 
error bound 

f d <(c,x)-(d,x-) + y-r=:J d . (31) 

Similarly as in the case of the lower bound, the upper bound uses the approximate value 
(c, x) and takes into account the violations of x wrt. to the cone % and the linear equations. 

The computation of the quantity (d,x~) can be avoided. Since x is an approximate optimal 
solution, x G % or x is close to % (provided the conic solver has computed reasonable 
approximations). Hence, for the supremum x + = sup{x, 0} the distance \\x — x + \\ is small. 
If we replace in Theorem 14.21 x by x + then the quantity |(c, x + ) — (c, x)\ is small, but 
x~ := < {x + , 0} yielding (d,x~) = and the upper bound 

fd< (c,x + ) + (y,r), (32) 

where \Ax + — b\ < r. In general it is not possible to compute the supremum x + exactly, but 
each close upper bound x> x + will suffice. 
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In the special case of linear programming we can take the exact supremum x + G R™ defined 
by (jHJ). Then x~ = 0, and we obtain the upper bound 

f<c T x + + y T r = J d . (33) 

The following short INTLAB program produces this upper bound: 

xtplus = max(0,xt) 

setround(-l) ; 

rn = abs(A*xtplus -b) ; 

setround(+l) ; 

rp = abs(A*xtplus -b) ; 

r = max(rn,rp) ; 

fu = c'*xtplus + yup J *r; 

setround(O) ; 

Until now we have assumed the existence of e-optimal solutions within some reasonable 
bounds. Now we mention briefly that in the case where appropriate primal or dual bound- 
edness qualifications are not known it is frequently possible to compute verified primal and 
dual feasible solutions which are close to optimality. These solutions can be used to compute 
verified reasonable error bounds for the optimal value. The basic algorithm consists of the 
following steps: 

(i) Perturb the original problem slightly such that the optimal solution of the perturbed 
problem is an interior feasible solution of the original problem. 

(ii) Solve the perturbed problem approximately. 

(iii) Use this approximation to compute an enclosure (i.e an appropriate interval) containing 
a feasible solution. 

(iv) Evaluate the objective function for the enclosure. 

Especially step (iii) is nontrivial since the existence of feasible solutions must be rigorously 
proved. Interval arithmetic provides several methods for computing enclosures of solutions 
for linear and nonlinear systems in the finite dimensional case, but also for infinite dimen- 
sional problems (certain types of ordinary and partial differential equations) enclosure meth- 
ods are known. However the bounds obtained in this way have two disadvantages. They 
are much more time-consuming than the previous ones, and they provide an upper bound 
of the primal optimal value only if the primal problem is well-posed, and a lower bound of 
the dual optimal value only if the dual problem is well-posed. 
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A detailed description of this algorithm can be found in the case of linear programming in 
[2], for smooth convex programming problems in [IB], and for semidefmite programming 
problems and linear matrix inequalities in [T6] . 



5 Second Order Cone Programming 

In SOCP the partial ordering is defined by the ice-cream cone @ yielding a finite dimensional 
vector lattice equipped with the following operations. 

Theorem 5.1 Let X C R n be defined by (TJj. Then for 16R" 

% j~ "JC fi 1 1 3j • 1 1 



and 



x^ = < 



X 







ry I , I rv> / ry> 
JU - I <Xj j i I Jb ■ 



I 2 b. 



x 



x 



if — \\x-\\ < x n < \\x : 



if x n — \\%: | 



if "^n — ||"^: 



\X. 



2 b. 



\x.\ 



if 



\X; < X n < \\X : 



(34) 



(35) 



Proof. First we prove ( 1341 . If x n > \\x-\\ then x E% which implies x + := sup{x, 0} = x. 
If x n — —\\ x :\\ then x G —DC and x + = 0. Finally, assume that — ||ar : || < x n < \\x.\\. Then 
X- 7^ 0, and a simple geometric argument shows that x + is the orthogonal projection of x 
onto the boundary of DC, that is 



x = a 



X- 



and = (x + — x) T x 



The latter condition describes the orthogonality. Since 



ax- 
a\\x-\ 

a 2 ||x.|| 2 



-a 

~Xrt. 



ax, 
a\\x- 1 



1 1 2 2 1 1 1 1 2 

Ot • I | Gl. ijC • I CC I • 



2a 2 \\x.\\ 2 — a(||x : || 2 + x n \\x.. 
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we get a = (\\x : \\ + x n )/2||x : || which proves ( l34"l) . 

Finally, ( 1351) follows from ( 134"1) since x~ = inf{x, 0} = — sup{— x, 0}. □ 

Due to rounding-off errors x~ may not be computed exactly. But for computing rigorous 
results we know from the previous section that it is sufficient to compute a lower bound 
xT < x~ by using directed rounding or interval arithmetic. The distinction of cases, however, 
must be implemented carefully when x n is almost equal to — or ||x : ||. A similar remark 
applies to the computation of an upper bound x + > x + . 

Let % be the Cartesian product of ice-cream cones %j C R n J for j = 1, . . . ,n. This is a 
convex, self-dual cone (see Alizadeh, Goldfarb [2]). The standard SOCP problem has the 
form 

n n 

minimize cjxj s.t. AjXj = b, Xj G %j for j = 1, . . . , n, (36) 

where Aj G R mxn i ; Cj,Xj G R nj and b G R m . If we merge these quantities 

A := (A\] . . . ; A n ), 

c := (ci;...;c n ), (37) 

x ■ (-^li • • • j X n )j 

then the standard SOCP problem has the form (|T4|) . and it follows that the dual problem 
( TTHj) can be written as 

minimize b T y s.t. — Ajy + Cj G %j for j = 1, . . . , n. (38) 

Here we have chosen the finite-dimensional spaces X := R™ where n = ^jUj and ^ := R m 
equipped with the Euclidean inner products. By x it j we denote the i-th component of the 
vector Xj, and x.j := (xij . . . ,x nj -ij) T . In this section 

X (S'li • • • ; 2-n) 

denotes a vector in % with x : j = and x n .j > for every j. Then 

x < x <^> ||a; :j -|| + x njJ < x nj)j for j = 1, . . . , n. (39) 

The computation of a rigorous lower bound for the optimal value of (1361) is a straightforward 
application of Theorem 14.11 

Corollary 5.1 Assume that PBQ holds for some x G %. Let y G R"\ and let 

dj := -Ajy + Cj for j = 1, . . . , n. (40) 
Suppose further that for j = 1, . . . , n there are lower bounds d~ < dj . Then: 
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(a) The primal optimal value is bounded from below by 

n 

f P >b T y + Y,d- 3j x nj ,-.= f_ p . (41) 
3=1 

(b) If d nji j > \\d;j\\ for j = 1, . . . ,n, then y is dual feasible and fd > / = b T y, and if 
moreover y is optimal then f d = f 

Proof. 

It follows from §M§ that 

f P > (y,b} + (d-,x) 

n 

= b T y+J2(ij,Xj) 

n 

= b T y+Y,d- jd x njd , 

where the last equation is fulfilled because x : j = 0. This finishes the proof of (a). If 
d nj ,j > \\d : ,j\\ for j — 1, . . . ,n then —Ajy + Cj G Xj, dj = 0, and y is dual feasible. Hence, 
(b) is proved. □ 

An upper bound for the optimal value of (1361) is an immediate application of Theorem 14.21 
Corollary 5.2 Assume that DBQ is fulfilled. Let x £ %, and suppose further that 

n 

I Yl ~ h \ - F ' 

i=i 

then: 

(a) The dual optimal value is bounded from above by 

n 

fd < ^2c j x j + y T f =: J d . 

3=1 

(b) If r — then x is primal feasible and f p < f d — Y^j=i c j^j> an d if moreover x is 
optimal then f p = f d . 

Proof. Since x £ %, it follows that x~ = < {x, 0}. Therefore the assertion is an 
immediate consequence of Theorem 14.21 □ 

SOCP solvers may compute an approximation x which is not contained in X, i.e. for some j 
the part Xj is not contained in the convex cone Xj. Then x is replaced by an upper bound of 
x + . As aforementioned, in floating point arithmetic formula (|34|) must be carefully evaluated 
using directed rounding such that the computed result is guaranteed to be in X. 
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6 Semidefinite Programming 

We examine the standard primal semidefinite programming problem 

minimize (C,X) s.t. (A u X) = b i: i = l,...,m, X e X, (42) 

where C, X and Ai are real symmetric sxs matrices, b G R m , X — SI, and ( , ) denotes the 
inner product ( |T3l) in the linear space of symmetric matrices. Using the svec operator (1121) 
such that 

c = svec(C), x = svec(X), a« = svec(v4j), (43) 
we can write problem (|42|) equivalently in the form 

minimize c T x s.t. Ax = b, x G X, (44) 

where A is the matrix with rows aj . Problem ( l44l) has the form (fT4l) . and therefore the dual 
problem ( [181) is 

maximize b T y s.t. - A T y + c G X, y G R' m . (45) 
The equivalent matrix notation is 

rn 

maximize b T y s.t. — s ^^yiAi + C E X. (46) 

2=1 

Corollary 6.1 Assume that the SDP satisfies: 

(i) Either the primal problem is infeasible, or 

(ii) there exists a nonnegative number x such that for every e > there exists a primal 
feasible solution X(e) <x-I and (C, X(e)} — f p < e, 

where I denotes the identity matrix. Let y G R m , and let 

in 

D = C-J2mA- (47) 

i=l 

Suppose further that dT < {\ min (D) , 0} , and that D has at most I negative eigenvalues. 
Then: 

(a) The primal optimal value is bounded from below by 

f P >b T y + l-d- -x=: f. (48) 



-v 



(b) If d =0 then y is dual feasible and fd> f — b y, and if moreover y is optimal then 

fd = L p - 



Proof. Let D have the eigenvalue decomposition D = Q T AQ, then (TTT1) implies D = 
Q T A~Q. Hence LT := d~ ■ Q T sign(-A~)Q < {.0,0}, where sign is +1,0 or -1 if the 
corresponding coefficient of the matrix is positive, zero or negative, respectively. Moreover, 
PBQ implies X := x ■ I > X (e) . Now from Theorem 14.11 (a) it follows that 

f P > (V,b) + (D- ,X) = b T y + 1 ■ d~ -x 

which proves (a). The part (b) is an immediate consequence of Theorem 14.11 (b). □ 

In order to control all rounding errors and to compute a verified lower bound / it is neces- 
sary to compute a rigorous lower bound of the smallest eigenvalue for a symmetric matrix. 
Interesting references for computing rigorous bounds of some or all eigenvalues are Floudas 
[5], Mayer [2T], Neumaier [22], and Rump [321 SO]- In VSDP we have computed the quantities 
I and d~ by using Weyl's Perturbation Theorem for symmetric matrices: For an approximate 
eigenvalue decomposition D = Q T AQ of D with eigenvalues A, on the diagonal of A, we use 
directed rounding or interval arithmetic for computing an error matrix E > | D — D \ . Then 
the Theorem of Weyl implies that 

\\{D)-\\ < \\E\\ 2 
for each eigenvalue \{D). Therefore, we obtain the bounds 

X- ||£||co < \{D) < A,+ H^Hoo (49) 
for all eigenvalues, yielding immediately the quantities d~ and I. The short INTLAB program 

[Qt.Lt] = eig (full (mid (D))) ; 
E = D - Qt * intval(Lt) * Qt ; ; 
r = abss (norm(E, inf ) ) ; 
lambda = midrad(diag(Lt) ,r) ; 

implements these bounds for the eigenvalues of a symmetric interval matrix D, where eig, 
abss and midrad denote the MATLAB and INTLAB routines for computing approximate 
eigenvalues, the absolute value, and the midpoint radius description of interval quantities, 
respectively. 

Corollary 6.2 Assume that DBQ is fulfilled. Let X e % = S+ and suppose further that 

\(A h X) -bi\ <fi fori = l,...,m. (50) 

Then: 
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(a) The dual optimal value is bounded from above by 

f d <(C,X} + y T f=:J d . (51) 

(b) Ifr — O then X is primal feasible and f p < f d — (C,X), and if moreover X is optimal 
then f p =J d . 

Proof. Because X~ = 0, this corollary is an immediate consequence of Theorem 14.21 □ 

In general, SDP-solvers, compute an approximation X which is not a positive semidefinite 
matrix, but there is a cluster of negative eigenvalues of X nearby zero. In order to enforce a 
positive semidefinite approximation, we compute a rigorous bound x < {A m i n (X),0}. Then 
it follows that X — xl G %, i.e. is positive semidefinite, and we can use this shifted matrix 
in the previous corollary. Then with directed rounding it is straightforward to compute f d . 



7 Block Structured Variables 

Frequently conic optimization problems have block structured variables, that is the variables 
are in the Cartesian product of different cones. More precisely, there are n real normed vector 
spaces Xi, . . . , X n , convex cones %\ C Xi, . . . , % n C X n , a real normed vector space V, and 
n continuous linear operators Aj : Xj — > V- Let X and % denote the Cartesian products of 
the spaces Xj and the cones X,-, respectively. The vectors x and c and the linear operator 
A are partitioned appropriately: 

x = (xi, . . . ; x n ) where Xj G Xj, 
c = (ci; . . . ; c n ) where Cj G X*, 
A = (A 1 ;...;A n ). 

Defining 

n n 

Ax:=^^AjXj and (c,x) :—'S^{cj,Xj), (52) 

it follows that A : X — > y is a continuous linear operator, and c G X*. The primal conic 
optimization problem with block structured variables has the form 

n n 

minimize ^j(cj, Xj) s.t. '^^AjXj = b, Xj G %j for j = 1, . . .,n, (53) 

and the dual problem is 

minimize (y, b) s.t. {-A\y\ -A* n y) + ( Cl ; . . . ; c n ) G %\ X . . . X X* n , y eW- (54) 
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The most important examples are semide finite- quadratic-linear programs. These are block 
structured problems where y = R' m , 

x — . . . , x ng , x 1 , . . . , x n? , x ; t a-, I^jj; 

and 

a;f , . . . , x s ns are symmetric matrices of various sizes, 

xf, . . . , x^ are vectors of various sizes, and 

x l is a vector. 

Using ( 1521) . ( 153|) . ( 1361) and ( )42|) we obtain the primal problem 

A\T<1 i (J\T r 



minimize ^ (c|, xf) + XI ( c t) x t + ( c ) x i 
j=i k=i 

n s rig 

s.t. E(4'.^ + E(4)M + (^ = ^ • = /,..., m 
i=i fc=i 



(56) 



x s 3 e%^ x q k e% q k , x l eX l Vj,k, 

where the matrices and vectors have appropriate dimensions, and Xj, % q k and % l are the 
convex cones of positive semidefinite matrices, the ice-cream cones, and the positive orthant, 
respectively. 

The dual problem has the form 

m 

maximize biUi 

i=l 

in 



s.t. V XjU; ■ <-~ e X* tor./ 1 n s 

i=l 

m 

-EAlyi + cleXl for* = l,...,n« 



(57) 



i=l 
in 



i=i 

Observe that the set of primal feasible solutions is the Cartesian product of semidefinite, 
quadratic and nonnegative orthant cones intersected with an affine subspace. It is possible 
that n s , n q or the length of x l is zero, which means that one or more of the three parts of 
the problem is absent. 

The following two corollaries provide finite lower and upper bounds of the optimal value for 
block structured problems. 

Corollary 7.1 Assume that PBQ holds for some x = (xi, . . . ;x n ) G X. Let y e ^, and 
assume that for j — 1, . . . , n 

dj := -A*y + Cj (58) 
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and dj <{dj,0}. Then: 

(a) The primal optimal value is bounded from below by 

n 

f P >(y,b) + J2(dj^)=:l p . (59) 

(b) If dj = for j = 1, . . .,n, then y is dual feasible and f d > f p — (y>b), and if moreover 
y is optimal then f d — f . 

Proof. This corollary follows immediately from Theorem 14. II by observing the linearity of 
the block structured variables. □ 

Corollary 7.2 Assume that DBQ holds for some y G y*. Let x = (xx, . . .] x n ) G % and 

suppose further that 

n 

\J2^j~b\<r, (60) 
i=i 

then 

(a) The dual optimal value is bounded from above by 

n 

TdZ^fazj) + f d . (61) 

(b) If r — then x is primal feasible and f p < f d , and if moreover x is optimal then 
fp = fd- 

Proof. The corollary is an immediate consequence of Theorem 14.21 □ 

Lower and upper bounds for the optimal value of semidefinite-quadratic-linear programs can 
be immediately obtained by inserting the preciding formulas for LP, SOCP and SDP. 



8 Certificates of Infeasibility 

A theorem of alternatives states that for two systems of equations or inequalities, one or 
the other system has a solution, but not both. A solution of one of the systems is called a 
certificate of infeasibility for the other which has no solution, since in principle this allows 
an easy check to prove infeasibility. Certificates of infeasibility are frequently computed 
by optimization algorithms if no feasible solutions of the primal or dual constraints exist. 
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Especially in the presence of equality constraints, certificates cannot be represented exactly in 
floating point arithmetic, and therefore approximate certificates can satisfy the constraints 
only within certain tolerances. This effect is amplified by rounding-off errors during the 
calculations for computing the approximate certificate. However, it turns out that in order 
to prove infeasibility by using floating-point arithmetic it is sufficient if an interval of small 
diameter can be given which guarantees to contain a certificate. We call such an interval a 
rigorous or verified certificate of infeasibility, and describe briefly how such certificates can 
be obtained for conic problems. We begin with two well known propositions and include the 
short proofs. 



Proposition 8.1 Suppose that y EY* satisfies 

A*yeX*, (y,b)<0, (62) 
then the system of primal constraints 

Ax = b, x G X (63) 

has no solution. 



Proof. If the system ( 163|) has a solution x, then < (A*y,x) = (y,Ax) = (y,b) contra- 
dicting our assumption (y, b) < 0. □ 

The linear functional y is called a certificate of primal infeasibility, and represents a dual 
unbounded ray. 

Proposition 8.2 Suppose that x G X satisfies 

Ax = 0, xeX, (c, x) < (64) 
then the system of dual constraints 

-A*y + ceX*, yeY* (65) 

has no solution. 



Proof. If the system ([65]) has a solution y G Y*, then < (—A*y + c, x) = ~{y,Ax) + 
(c, x) = (c, x) < contradicting our assumption. Hence, system (16"5"]) has no solution. □ 

The vector x is called a certificate of dual infeasibility and represents a primal unbounded 
ray. 

Many conic solvers expose infeasibility by computing approximate unbounded rays. Given 
an approximate primal unbounded ray x G X, dual infeasibility is proved if the equation 
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and sign conditions (164j) can be checked rigorously on a computer. The underdetermined 
equation Ax = is in general not exactly satisfied for floating-point certificates x. To 
obtain a rigorous certificate we proceed as follows: Let (3 be an approximation of (c, x), 
and assume that (3 < 0. Otherwise, for nonnegative (3 the sign condition in f[6"4"j) would be 
even not satisfied for the approximate primal unbounded ray x, and this indicates that x is 
not suitable. Then we compute an interval of small diameter x, also called enclosure, for a 
solution of the underdetermined linear system 

Ax = and (c, x) = (3 < 0, (66) 

which is close to x. If x C then there exists an x e x which fulfills the condition (1641) 
yielding a rigorous certificate of dual infeasibility. This check depends on the special problem 
and requires further information about the operator A and the cone %. In the three cases LP, 
SOCP, and SDP for the finite-dimensional linear system (I6"6"j) methods of interval arithmetic 
can be used for computing an appropriate enclosure x. For a detailed description of such 
an algorithm see [H]. The condition x = [x,x] C % can be verfied for LP by checking the 
equivalent condition 

x > 0, (67) 
for SOCP we check the equivalent condition 

X-n > ( 68 ) 

and for SDP we check 

A min (X) > 0, (69) 

where x = svec (X) . 

In the case of an approximate dual improving ray y, primal infeasibility can be rigorously 
proved on a computer if the condition ( 1621) can be verified; that is, if the sign conditions 

(y, b) < and (A*y, x) > for all x E % (70) 

can be checked reliably. As before, this check depends on the special problem and requires 
further information about the operator A* and the cone %*. 

It follows immediately that for LP (170|) is equivalent to 

b T y<0 and A T y > 0, (71) 
for SOCP we obtain the equivalent condition 

b T y<0 and (Ajy) nj > \\(Ajy) : \\ for j = 1, . . . , n, (72) 
and for SDP we get 

m 

b T y<0 and A min (^^) > 0. (73) 
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All three conditions can be checked rigorously by using directed rounding or interval arith- 
metic. The vector y provides a rigorous certificate which can be viewed as a degenerate 
interval with zero diameter. 

9 Combinatorial Optimization 

Linear and semidefinite programs play a very useful role in global and combinatorial op- 
timization (Wolkowicz [4"7]). Several methods (for example lift-and-project methods) are 
known for constructing linear or semidefinite relaxations, which are used in branch-bound- 
and-cut algorithms to eliminate regions which do not contain global minimizers. Neumaier 
and Shcherbina [31] have pointed out that backward error analysis has no relevance for com- 
binatorial programs, since slightly perturbed coefficients no longer produce problems of the 
same class. There, one can also find an innocent-looking linear integer problem for which 
the commercial high quality solver CPLEX [12] and several other state-of-the-art solvers 
failed. The reason is that the relaxations are not solved with sufficient accuracy and global 
minimizers are truncated. Hence, in order to obtain safe results, it is important to have 
reliable, good and cheaply computable lower bounds of the optimal value for relaxations. 

Various problems like Max-Cut, Partitioning, Coloring and many others can be formulated 
as linear integer problems where the vector of decision variables x G {—1,1}". Tight semidef- 
inite relaxations are obtained by lifting the vector x into the space of semidefinite matrices 
by the operation 

X = xx T . (74) 

It follows immediately that 

X y 0, diag(X) = e, and rank(X) = 1, (75) 

where e is the vector of ones. Dropping the condition rank(X) = 1 we obtain a semidefinite 
relaxation. Laurent and Poljak [20] have shown that for this type of relaxations — 1 < Xij < 
1, and if X^ G { — 1, 1} then X = xx T where x G { — 1, l} n . This property establishes the 
tightness of these relaxations. Moreover, it follows that the primal boundedness qualification 
is fulfilled in the way that an optimal solution exists with A max (X) < n, and thus rigorous 
lower bounds for the optimal value can be computed. 

Sometimes, these tight relaxations are in addition ill-posed. As an example we consider 
Graph Partitioning Problems. These are known to be NP-hard, and finding an optimal 
solution is difficult. Graph Partitioning has many applications among those is VLSI design. 
Here, we investigate semidefinite relaxations for the special case of Equicut Problems, which 
have turned out to deliver tight lower bounds (see also Gruber and Rendl [H]). The general 
case of Graph Partitioning Problems can be treated similarly. 
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Given an edge-weighted graph G with an even number of vertices, the problem is to find a 
partitioning of the vertices into two sets of equal cardinality which minimizes the weight of 
the edges joining the two sets. The algebraic formulation is obtained by representing the 
partitioning as an integer vector x G { — l,l} n satisfying the parity condition ^ /i Xi = 0. 
Then the Equicut Problem is equivalent to 

1 — x x n 
min / ^ ajj - J subject to x E { — 1, 1}™, ^J^j = 0, 

i<j i=l 

where A = (a^) is the symmetric matrix of edge weights. This follows immediately, since 
1 — XiXj = iff the vertices i and j are in the same set. The objective can be written as 

- Qjj(l — XiXj) = -x T (Diag(Ae) — A)x = -x T Lx, 

Zi 4r 4r 

i<j 

where L := Diag(Ae) — A is the Laplace matrix of G. Using x T Lx = trace(L(xx T )) and 
X = xx T , it can be shown that this problem is equivalent to 

f p = min-(L, X) subject to diag(X) = e, e T Xe = 0, X >z 0, rank(X) = 1. 

Since X >z and e T Xe = implies X to be singular, the problem is ill-posed, and for 
arbitrarily small perturbations of the right hand side the problem becomes infeasible. By 
definition, the Equicut Problem has a finite optimal value f p , and a rigorous upper bound 
of f p is simply obtained by evaluating the objective function for a given partitioning integer 
vector x. Hence, it is left over to compute a rigorous lower bound. At first, the nonlinear 
rank one constraint is left out yielding an ill-posed semidefinite relaxation, where the Slater 
condition does not hold. The related constraints diag(X) = e and e T Xe = can be written 
as 

(Ai, X) = bi, bi = 1, Ai = Ei for i = 1, . . . , n, and A n+1 = ee T , b n+1 = 0. 

where Ei is the n x n matrix with a one on the ith diagonal position and zeros otherwise. 
Hence, the dual semidefinite problem has the form 



n 

max 



- 1 

s.t. diag(j/i.n) + y n +i(ee T ) ^ -L, ye R n+1 . 

i=l 

The constraints diag(X) = e, X >z imply PBQ with finite upper bounds A max (X) < x = n. 
Corollary 16.11 yields 

Corollary 9.1 Let y e R n+1 , and assume that the matrix 

D = -L - Diag(y! . n ) - y n+1 (ee T ) 
has at most I negative eigenvalues, and let d < X m i n (D). Then 

n 

fp >^Vi + l ■ n ■ dT -■ f_p 

i=l 
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n 


L 


M/p> Id) 


Kfd, l p ) 


t 


tiow 


tc 


100 


-3.58065e+003 


-7.117e-008 


3.843e-011 


4.2 


0.5 





200 


-1.04285e+004 


-7.018e-008 


9.621e-010 


7.9 


0.2 





300 


-1.90966e+004 


-2.573e-008 


8.918e-009 


21.1 


0.9 





400 


-3.01393e+004 


-1.633e-008 


3.008e-008 


39.0 


2.0 





500 


-4.22850e+004 


1.431e-008 


2.584e-008 


67.5 


3.7 





600 


-5.57876e+004 


5.418e-009 


1.829e-008 


124.7 


6.0 


-5 



Table 1: Results for Graph Partitioning 

In Table Q] some numerical results for problems given by Gruber and Rendl pT] are displayed. 
The number of nodes is denoted by n. For this suite of ill-posed problems with up to 601 
constraints and 180300 variables the semidefinite programming solver SDPT3, version 3.02 
[4*6] has computed approximations of the dual optimal value fd, which are close to the 
approximate primal optimal value f p ; see the column fi(f p ,f d ). Here, the relative accuracy 
of two real numbers a and b is measured by the quantity 

u 

n(a, b) :- 



max{1.0, (\a\ + \b\)/2}' 



The negative signs in this column show that weak duality is violated for the computed 
approximations in four cases. SDPT3 gave tc = (normal termination) for the first five ill- 
posed examples. Only in the last case n = 600 the warning tc = —5 (that means : Progress 
too slow) was returned. We have computed the lower bound / by using Corollary 19.11 The 



— p 



small quantities /i(/d, / ) show that the overestimation of the rigorous lower bound / can 
be neglected. In Table [T] the times for computing the approximations with SDPT3, and for 
computing / by using Corollary 19.11 are denoted by t and ti ow , respectively. It follows that 
the additional time ti ow for computing the rigorous bound / is small compared to the time t 
needed for the approximations. This is of the same tenor as the quotation of Turing at first. 



10 Numerical Results 



In this section we describe briefly our numerical experience. Lurupa [IT] is a C++ imple- 
mentation of the presented rigorous bounds for the special case of linear programming. In 
the following we give a short summary of numerical results for the NETLIB suite of linear 
programming problems [26]. For details refer to p~9]. The NETLIB LP-suite is a well-known 
collection of difficult to solve problems with up to 15695 variables and 16675 constraints. 
They originate from various applications, for example forestry, flap settings on aircraft, and 
staff scheduling. We chose the set of problems for which Ordonez and Freund [33] have 
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computed condition numbers. There it is stated that 71% of the problems have an infinite 
condition number. As Fourer and Gay [7] observed, preprocessing can change the state of an 
LP from feasible to infeasible and vice versa, and therefore preprocessing was not applied. 

Roughly speaking, a finite lower bound (upper bound) of the optimal value can be computed 
iff the distance to dual infeasibility (primal infeasibility) is greater than zero. For 76 problems 
a finite lower bound could be computed with a median accuracy of med(/x(/, /_)) = 2.2 ■ 10~ 8 
(/ is the approximate optimal value) and a median time ratio of med(ti ow /t) = 0.5. For 
35 problems Lurupa has computed a finite upper bound with med(fi(f,f d )) = 8.0 ■ 10~ 9 
and med(t up /t) = 5.3. For 32 well-posed problems finite rigorous lower and upper bounds 
could be computed with med(/z(/ p , /,)) = 5.6 ■ 10" 8 . Only for two problems with finite 
condition number (SCSD8 and SCTAP1) an upper bound of the optimal value could not 
be computed. Taking into account the approxmimate solver's default stopping tolerance of 
10 -9 , the guaranteed accuracy computed with Lurupa for the NETLIB LP suite is reasonable. 
The upper bound is more expensive, since linear systems have to be solved rigorously, and 
sometimes perturbed problems have . 

For the SDPLIB benchmark problems of Borchers [I] we have comuted with VSDP [TS], a 
MATLAB software package for verified semidefinite programming, rigorous bounds . For de- 
tails see [15] and [16]. Freund, Ordonez and Toh [8] have solved 85 problems with SDPT3 out 
of the 92 problems of the SDPLIB. They have omitted the four infeasible problems and three 
very large problems where SDPT3 produced out of memory. In their paper interior-point 
iteration counts with respect to different measures for semidefinite programming problems 
are investigated, and it is pointed out that 32 are ill-posed. VSDP could compute (by using 
SDPT3 as approximate solver) for all 85 problems a rigorous lower bound of the optimal 
value and could verify the existence of strictly dual feasible solutions, which proves that all 
problems have a zero duality gap. A finite rigorous upper bound could be computed for all 
well-posed problems with one exception; this is hinf 2. For all 32 ill-posed problems VSDP 
has computed the upper bound f d = +oo, which reflects exactly that the distance to the 
next primal infeasible problem is zero as well as the infinite condition number. 

For the 85 test problems (not counting the 4 infeasible ones) SDPT3 (with default values) 
has computed good approximations and gave 7 warnings, and 2 warnings are given for well- 
posed problems. Hence, no warnings are given for 27 ill-posed problems with zero distance 
to primal infeasibility. In other words, there is no correlation between warnings and the 
difficulty of the problem. At least for this test set our rigorous bounds reflect the difficulty 
of the problems much better, and they provide safety, especially in the case where algorithms 
subsequently call other algorithms, as is done for example in branch-and-bound methods. 

Some major characteristics of the numerical results of VSDP for the well-posed SDPLIB- 
problems are as follows: The median of the time ratio for computing the rigorous lower 
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(upper) bound and the approximation is 0.085, (1.99), respectively. The median of the 
guaranteed accuracy for the problems with finite condition number is 7.01 • 10~ 7 . We have 
used the median here because there are some outliers. 

One of the largest problems which could be rigorously solved by VSDP is thetaG51 where 
the number of constraints is m = 6910, and the dimension of the primal symmetric matrix 
X is s = 1001 (implying 501501 variables). For this problem SDPT3 gave the message out 
of memory, and we used SDPA [9] as approximate solver. The rigorous lower and upper 
bounds computed by VSDP are = -3.4900 • 10 2 , J d = -3.4406 • 10 2 , respectively. This is 
an outlier because the guaranteed relative accuracy is only 0.014, which may be sufficient in 
several applications, but is insufficient from a numerical point of view. However, existence 
of optimal solutions and strong duality is proved. The times in seconds for computing 
the approximations, the lower and the upper bound of the optimal value are t= 3687.95, 
tiow — 45.17, and t up = 6592.52, respectively. 

For further numerical results and applications concerning ill-posed problems and the problem 
of computing the ground state energy of atomic and molecular systems by using a variational 
approach (see for example Fukuda et al. [TU] and Nakata et al. [23]) refer to VSDP [To] . 

To our knowledge no other software packages compute rigorous results for semidefinite pro- 
grams. There are several packages that compute verified results for optimization problems 
where the objective and the constraints are defined by smooth algebraic expressions. Elab- 
orate comparisons with some of these packages in the case of linear programming problems 
can be found in the forthcoming paper of Keil [TS] . 

11 Conclusions 

The computation of rigorous error bounds for conic optimization problems can be viewed as 
a carefully postprocessing tool that uses only approximate solutions computed by any conic 
solver. The bounds are developed in the framework of functional analysis. Error bounds for 
special conic problems can be derived easily. 

Several numerical results demonstrate that rigorous error bounds can be reasonably easily 
computed even for problems of large size and for ill-conditioned problems, in most cases with 
a range which is not much larger than necessary. 
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